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Abstract 

We present a combined theoretical approach to study the nonequihbrium transport properties 
of nanoscale systems coupled to metallic electrodes and exhibiting strong electron-phonon inter- 
actions. We use the Keldysh Green function formalism to generalize beyond linear theory in the 
applied voltage an equation of motion method and an interpolative self-energy approximation pre- 
viously developed in equilibrium. We analyze the specific characteristics of inelastic transport 
appearing in the intensity versus voltage curves and in the conductance, providing qualitative cri- 
teria for the sign of the step-like features in the conductance. Excellent overall agreement between 
both approaches is found for a wide range of parameters. 

PACS numbers: 73.63.-b, 71.38.-k, 73.63.Kv 



I. INTRODUCTION 



Advances in the field of molecular electronics and nano-objects pQ have motivated an in- 
creasing interest in electron-phonon interaction [2j. Experiments give evidence that electron- 
vibrational coupling within the molecule play an important role in its charge transport prop- 
erties. This was first found by Park et al. on Cqq [3J. Also, the excitation spectra show 
features that could be ascribed to sidebands formed by the presence of strong electron- 
phonon interactions [U E] . 

From the theoretical point of view, the problem of the interaction of a localized level with 
a field of bosons can be traced back to the small polaron model of Holstein [6J . Today the so- 
called Anderson- Holstein Hamiltonian is the simplest and more commonly used Hamiltonian 
to study the electronic transport through molecular systems. This Hamiltonian has not an 
exact solution except for a few special cases in equilibrium. Therefore it is desirable to 
develop different theoretical approaches which would allow to calculate and predict robust 
behaviors for physical magnitudes directly comparable to out of equilibrium experiments. 
With this aim we use a Keldysh Green function formalism [7J to generalize two theoretical 
approaches previously developed by us in equilibrium, the equation of motion (EOM) method 
[8] and the interpolative self-energy approximation (ISA) [9] , to deal with situations in which 
many phonons can be absorbed/emitted by the molecular system when a bias voltage is 
applied between the electrodes. This is clearly a nonequilibrium situation which cannot be 
described by extensions of equilibrium theories to small voltages, if the voltage exceeds the 
phonon frequency. 

Previously, the problem of the electronic transport through molecular junctions or quan- 
tum dots has been approached in different ways depending on the different regimes deter- 
mined by the parameters: the temperature T, the coupling of the localized level to the leads 
characterized by the level width P, the coupling of the localized level to phonons A, and 
the phonon frequency ojq. The semi-classical regime, defined by T >> P, can be described 
from a master equations point of view [TUl E]. In the quantum regime T << P, the ra- 
tio A/P distinguishes between the weak and strong coupling regimes. The weak coupling 
regime, A/P << 1, can be approached by a variety of methods with the common charac- 
teristic of being perturbative in A/P such as the Born approximation or the self-consistent 
Born approximation [T2HTU [161 126], perturbative renormalization theory [IT] or diagram- 
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matic techniques [HI [181 IH]- this respect we should mention that the two approaches 
introduced in the present work recover this hmit. The quantum, strong couphng regime, 
A/r >> 1, for which perturbation theory breaks down is much more difficult to analyze and 
the decoupling of electronic and vibronic degrees of freedom has been a usual approximation 
[T6l I20ti22] . This work concentrates in this limit approaching the problem from two very 
different starting points. 

A remarkable experimental result is the observed step-like feature in the differential con- 
ductance at bias voltages equal to the phonon energy that can be either upwards or down- 
wards [SHSl I23ti26] . This fact has attracted a considerable theoretical interest [21 [13] lately. 
In the limit uq << F, a symmetric contact and small A the behavior turns out to depend 
only on the transmission r of the junction with the step upwards (downwards) for r < 1/2 
(r > 1/2) [2S1 [2Z|- This result seems to offer a rough rule of thumb for predicting the 
observed step sign. Outside this limiting situation this feature on the conductance will de- 
pend in a more complicated form on the system parameters [281 ES] • The same issue will be 
addressed in this work in the strong coupling regime. 

In order to introduce the method we will consider in this paper the spinless version of 
the Anderson- Holstein model [30ti32] . In section I we introduce the nonequilibrium Green 
functions formalism used for the calculation of the transport properties of the system. In 
sections II and III we present the out of equilibrium extensions of the EOM method and the 
ISA respectively. Section IV is devoted to the analysis of the intensity versus voltage curves 
and the conductance as obtained by both approximations. The remarkable overall agreement 
found between such different theoretical approaches in the out of equilibrium situation for a 
wide range of parameters, gives confidence in our results. We find the I-V curves to increase 
stepwise when a new inelastic channel emitting n phonons opens. The conductance reveals 
more interesting features of the emission processes. While its main peak, obtained at low 
voltages, is almost identical to the main resonance appearing in the equilibrium density of 
states, the phonon side-bands show specific behavior associated to inelastic transport which 
therefore cannot be obtained by any extension of equilibrium calculations to finite voltages. 
The origin of such features is analyzed. Finally, our conclusions are presented in section V. 
Atomic units e = h = m = 1 are used throughout this work except otherwise stated. 
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II. GENERAL NONEQUILIBRIUM FORMALISM 



We consider the spinless Anderson-Holstein Hamiltonian describing a single non-degenerate 
electronic level, eo, coupled linearly to a local phonon mode of frequency Uq and to electronic 
reservoirs, 

H = eocjco + (yk,uc\ ,,CQ + c.c. j + ^{ek,v + ^-y)c\^^Ck,y + Wo^"^^ + + 6)no (1) 

k,v k,v 

where is the phonon energy, A the electron-phonon coupling constant, e^,^ with u = L,R 
denotes the single particle energies of the left and right electrodes, fiL — fJ'R = being the 
applied bias and Vk^u the coupling between the localized level and the reservoir states. 

The electronic transport properties through this system can be conveniently calculated 
using the nonequilibrium Green function formalism or Keldysh method |7]. For a stationary 
situation the retarded and the nonequilibrium distribution Green functions G'^~ and 
G ^ are defined as follows: 

G:^{co) = -t j e{t - t') < c]{t)c,{t') + c,{t')c]{t) > e'^^'-''^d{t - t') 
Gtf{uj) =^J< c]{tMt') > e*'^(*-*')rf(t - t') (2) 
Gr^+{io) = -^J< c.{t')c]{t)) > e*-(*-*')rf(t - t') 

The frequency dependent Keldysh Green functions can be obtained from the correspond- 
ing Dyson equations which in matrix form read: 



G+- = + g+"S'^G^ + g'"S'"G+- - g^'S+'G^ (3) 
G+- = g+- + G+^SV" + G''S'"g+- - G'S+^g"^ 
G+- = (I + G'"S'")g+-(I + S^G"^) - G^'S+^G^ 

where I is the unit matrix and g are the Green functions of the uncoupled system {Vk^u = 0) 
and with similar equations for the G"^ functions. The crucial point within this formalism 
consists in finding a reasonable approximation for the self-energies S*" and The current 
intensity between the reservoir u and the quantum level can be written in terms of the G^ 
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Green functions as: 

= ^ E Vu^. I du[Gl-,{u) - GtM] (4) 
k 

where the subindex labels the dot level. 

Using Eqs.(|4]) it is possible to write the current density in terms of the dot level Green 
functions. In particular, Eqs. (3) lead to |33] : 

k k 

where the Green function Gq^ is calculated from the corresponding Dyson equation: 

^00 = X] \ ^k,u\'^\Goo\'^g^f^^^ - IGooP^oo (6) 

k,u 

with a similar equation for Gq^^ . All the expressions can be simplified by making the usual 
wide-band approximation [31]: 

El^'c1'C.H = 2^r./-M (7) 

k 

where fui^^) are the Fermi distribution functions of the electrodes and T^, are taken as 
constants. 

From Eqs.(|5| and ^ the current can be written as a sum of an elastic and an inelastic 
contribution, 1^ = ijf''^ + lu"^^ as: 



= -^r. J rfa;|GSoMrPoo^M/.M + SoVM(i - /.M)] (8) 

Due to current conservation, II = —Ir, an equivalent expression can be obtained by 
means of the identity / = {TrIl — TlIr)/T with T = Tl + leading from Eqs. ([s]) to the 
well known expression [31] : 

/ = I duImGUu)[h{u^) - Uuj)] (9) 

From Eq.([9]), the differential conductance is obtained as G = dl/dV. In the linear 
response regime V ^ 0, Gqo{u) can be evaluated in equilibrium, G^oiA — Gqq'^{u), and the 
conductance can be expressed in terms of ImGQQ'^{u = (j^l) and ImG'^Q'^iuj = (Ir). However, 
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this is not in general a good approximation for V > uq and a full calculation of Gqq^u) has 
to be performed. 

On the other hand, the level occupation can be obtained from the non equilibrium spectral 
density functions as: 

1 



2m 
1 
27ri 



duG^Q (w) = < n > 



dcoFooico) = 2<n> -1 (10) 



where Fqq = Gqq + Gqq 



III. EQUATION OF MOTION METHOD FOR THE ANDERSON-HOLSTEIN 
HAMILTONIAN OUT OF EQUILIBRIUM 

In the quantum strong coupling regime we are interested in, it is convenient to apply 
to Hamiltonian Eq.(l) a standard canonical transformation H = SHS^^ with S given by 

[SHIES] 

S = exp[—{b^ - b)ho] (11) 
which transforms electronic and bosonic operators as 



Co = Coexp[- — {b^ - b)] 



b 



b no 



(12) 



Note that Eqs.(12) imply that the number operators for electrons in the level and in the 



leads remain unchanged. Then, the transformed Hamiltonian reads: 



H = eouo + ^ ek,„hk,^ + ^ Vk,^{cl ,^Co + clck,^) + uj^b^b (13) 

k,v k,u 

with io = eo — X^/uq representing the renormalization of the energy level due to its coupling 
with the local phonon. 

The nonequilibrium Green's functions will also be written in terms of the tilde-operators. 
In the EOM procedure, we will obtain Green's functions for other operators 0{t) different 
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from co(t) at time t, which are defined in a way similar to Eqs.(|3]) but with a more convenient 
notation. Also, instead of the functions G^^ and G ^ it is more convenient to use here their 
sum F. Then we write 



G^iO; t, t') = i9{t' -t)< d{t)cl{t') + cl{t')d{t) >fj 
F{d;t,t') = I < 5l{t')d{t) - d{t)i{t') >^ (14) 

From now on, the symbol < .... >fj means that the average should be taken with respect 
to the transformed Hamiltonian H. 

The EOM method for solving the Anderson-Holstein Hamiltonian was already introduced 
in [8] . Briefiy, starting with G"(co; t, t') from Eq.([l4)) and applying the equation of motion, a 
hierarchy of new Green's functions G°'{h'^^CQV ]t,t') is generated. To obtain a closed system, 
at a given step of the procedure we contract pairs of operators Ck^u and c\, ^ where possible 
as 

c\iyCk,u ~ Sk,k'SuM' < nk,u > (15) 

In this equation < 71^^ > is the Fermi-Dirac distribution function of the i/-electrode. Due to 
the fact that the non equilibrium problem we are interested in is much more involved than 
the equilibrium one addressed in |8j, we will restrict the method to the order 0{V^j^). Then, 
the following system of linear equations has to be solved: 



- eo - r{uij))G^{b^'%bi;uj) =< dlb^%b^ + h'^'c.Vcl >^ + 



Wo 



i-i 



XG^ib^'doV-^^; to) + XG^ib^'^^CoV; u) + 



1=0 



m=l 



m j \ I 



1=0 V^O/ rn=l 



i—m I II I 



m j \ I 



mj ) 



(16) 



n 
I 



with 1 I = lUnLiv ■ We have defined ujij = uj + {i — j)uo and the advanced self-energies 



T(uij) = y 

(^ij - ^k,u - if] 

k,u •' ' 

V^T^2 <nk,u> 



kM •' ' 



T] being an infinitesimal. In the wide-band limit to be used in this work T(ujij) = i(rL + Tji). 

We should point out that this procedure does not decouple electrons and phonons as it 
has been frequent in the literature. Rather, quantum coherence is preserved in all of the 
Green functions G"^(6^*Co6-'; uj) which involve emission of i and absorption of j phonons. On 
the other hand, since we have decoupled the Green functions involving the localized level and 
the electrodes to the order 0(V^^^), the procedure is somehow perturbative in Vk^u- However 
it becomes exact not only in the limit Vk^,y — ?■ but for A —t- and finite Vk^u as well. 



In Reference fS] we argued that all the expectation values of the type < b'^^clck^ub"^ > 



H 



appearing in Eq.(16) could be neglected. This is not in general the case when an electric cur- 
rent circulates through the localized level because these expectation values just describe the 
transit of an electron from the electrode to the level with absorption of m and emission of n 
phonons, which is the process we are analyzing. Therefore, they have to be calculated consis- 
tently with the appropriate non equilibrium Green's functions F's as we will explain below. 
With respect to < 6"''"c|)Co6™ >/^, these expectation values describe fiuctuations in level oc- 
cupancy when phonons are absorbed and emitted and should also be calculated consistently 
with the appropriate F's functions. However, we have checked that the approximation 



< fet-gtgofo- > 8^o5n^ < 4co > (18) 

is still a good approximation out of equilibrium at zero temperature. 

The calculation of the F's Green's functions follows the same lines even though it is more 
involved. Starting from F{co]t,t') and applying the equation of motion, we obtain new 
Green's functions, which are calculated from their equations of motion. A typical equation 
being 
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dF{b^%b^;t,t') 



dt 



-i{~eo - {t - j)uo)Fib^%bi;t,t') - 

-YclcoCk..^; t, t') + Fib^%~clckAb + -Y; t, f)] 



(19) 



In the next step, the F's functions appearing in the forth term of Eq.(19) are calculated from 



their EOM and approximated by the contraction of operators indicated in Eq.(15), yielding: 



^F(6t"gtgoCfc,,6'";t,t') = -z(efc,, - (n - m)coo)F{b^^5lcoCk,ub^;t,t') + 

iVk,u < nfc,. > F((6t + -fc^b"'- 1, t') (20) 



and 



dt 



F{b^''do5lck,,b"';t,t') = -i{ek,. - (n - m)uo)F{b^''do5lck,ub"';t,t') + 



iVk,u < 1 - Uk,^ > F(6"^"co(& 



(21) 



Eqs.(20) and (21 ) are now integrated in time from an initial time t = to where the system 



starts to evolve, with the initial conditions 



F{b^^5lcoCk,,b"';to,t') =< 2nu,u - 1 > G\b'^'^clc^cu,,b'^M') 



(22) 



and 



F(fc^"£o4cfc,.&'";to,t') =< 2n,,, - 1 > G\b^-^coclcu,r,b'^;h,t') 



(23) 



These equations come from the general definitions of Eq.( 14) by taking into account that. 



initially, the localized level and the leads were non-interacting independent systems. Since 
the EOM for the advanced Green's functions were previously derived, they are integrated 
backwards in time t, from its final value t' to its initial value to- Then we obtain 



A 



rt' X 
iVk,u<nk,,> [ dTF{{b^ + -rdob'--T,t')e 

Jtn ^0 



+ 



(24) 



and 



A 



-t' 

^m. ^ ^'■j^-ii<ik,i^-{n-m)uio){t-r) _j_ 

'to 



iVkM < 1 - rife, 



J to 

iVk,u < 1 - > / dTF{h^''cQ{h - 

Jtn 



A . 

UJo' 

A , 



Wo 



m. ^ j.l^^-i(ek,u-(n-m)u)o){t-T) 



(25) 



When to ^ — oo, the Fourier transform of Eqs.(24) and (25) can be readily obtained after 



taking into account that the integrals appearing in these equations can be written as the 



convolution product of two functions. Eq.(19) is also Fourier transformed yielding the final 



expression that allows us to obtain the F's Green's functions from: 



LiJij — Co 



r*(u;,,))F(6t*go&'';t^) = VL{uJ^J)G''{h^%V■uJ) + 



1=0 



A A A 

1^0 / ^ ^0 



j-l 



A 



i27r V f I (—\ V Vk.u < 2nk,^ - 1 >< (6^ - — )*4cfc,i.&' >h ^i^u ' ^k,u) + 
1=0 \U ^^o/ ^ Wo 



AF(6t^go6i+i; w) + AF(6t(*+i)go6J; co) + 



j-l 



i-1 

z=o 



A 



i-l 



[G%b^%b';cu)J2i-^) 

m=l 

F{b^%b';oo)j2i-^) 

m=l 
i 

[G^{b^'coV;u)Y,i-^) 

m=l 

i 

F{b^%b^;u)Y,{-iy- 



j-m 



m 



m \ I 



m 



m \ I 



m 



m \ I 



m 



m=l 



m \ I 



r^'>{u^,)] (26) 



10 



where we have defined the following self-energies: 



^{uij) = i2Tx ^ < 2nk,u - 1 > S{ujij - ek,u) 
Vt''''\ujij) = i27c ^ V^^ < 2nk,u - 1 >< rik^u > S{ujij - et.u) 
n^''\io^j) = i27T J2 < 2^fc.^ -1X1- Uk,^ > 6{uj,j - ek,u) (27) 



k,i/ 



The linear sets of Eqs. (16) and (26) are coupled trough the different expectation values 
appearing in these equations, which have to be calculated self-consistently. To do so, notice 
that from the definitions of F{b^"-clcoCk,ub'^]t,t') and F{b^"'Cob"']t,t') and for t = t' ^ +oo 
one has the identities 



^F(6t"gt5oCfc,,6'";u;) = F{b'^'-clUk,ub''\t\t') = i < 6t"4cfc,,(6 + >^ (28) 



and 



fj, , A A 

— F(6t"go6™; oo) = F(6t"ao6-; t', t') = i< 6^"co4(fe + + (b^ - — )"4co&" >h 

-oo 27r LUq OJq 

(29) 

respectively. By making use of these relations we can obtain the required expectation values 



from the EOM of the F's Green's functions. Once the system of Eqs. (16) and (26) are 
solved, the current I is calculated from Eq.([9]). An important point is related to current 
conservation, 1^ = —Ir, which is not automatically satisfied for a given approximation (see 
[571 ES])- We have numerically checked that the EOM method fulfills current conservation 
within the accuracy of the calculation, in the range of parameters investigated in the present 
work. 



IV. INTERPOLATIVE SOLUTION FOR THE ANDERSON-HOLSTEIN HAMIL- 
TONIAN OUT OF EQUILIBRIUM 

In this section we will introduce an interpolative approach for the calculation of the self- 
energy out of equilibrium. This approach is a generalization of a previous one developed 
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for an equilibrium situation [9]. It has also been successfully applied to a purely electronic 
problem like a quantum dot out of equilibrium |38j . An interpolative approach is possible due 
to a property of the self-energy which exhibits the same mathematical form when expanded in 
the interaction parameter [391 - HT] (which in Hamiltonian Q is the electron-phonon coupling 
A) both in the atomic [Vk^u 0) an in the perturbative (A — ?■ 0) limit. 

We briefly summarize the interpolative approach in an equilibrium situation. In the 
\4,, ^ hmit Eq.Q can be exactly diagonalized by means of a canonical transformation 
yielding for the level retarded Green function: 



a2 \2m 



Gro'iu) = e ^ > r + 30 

(jOqUiI \uj — eo — mujQ + 11] u — eo + mujQ + irj J 



where eo = cq — A^/cuq and (n) is the level occupation. From Eq. (30) we can calculate the 
expression for the level self-energy by means of the corresponding Dyson equation Sqq*^ = 
u — en — Gqo*'' ^ where en = — 2(A^/a;o)(?T.) is the energy level corrected by the Hartree 
contribution. In the limit of small electron-phonon coupling X/uq << 1 and to order A^, the 
atomic self-energy tends to: 

E£V)^A'r + J (31) 



On the other hand, the retarded self-energy of this model can be calculated up to A^ from 
the appropriate diagrams using perturbation theory [9J. In addition to a constant Hartree 
contribution this self-energy has the form: 

where p^^\uj) = F/ [(w — Cg//)^ + F^] /vr is the level density of states of the one-electron 
unperturbed problem, eg// being an effective level position which can be used for achieving 
charge consistency between the one-electron and the interacting cases (see for details). 
In the limit F — )• the above expression tends to 

22\q'{u) A = F(uj) (33) 

The interpolative self-energy is then calculated by means of the following ansatz: 
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2^00 ('^ 



)] } (34) 



where is the inverse function defined by Eq. (33). This ansatz recovers both the atomic 
hmit r/A — 7- and the opposite hmit where perturbation theory is vahd A/F — )■ and is 
in excellent agreement with NRG calculations and exact finite system diagonalizations in 
parameter space |9]. 

This ansatz can be generalized for a nonequilibrium stationary situation like the one 
addressed in the present work (previous theoretical approaches have been restricted so far 
to the case of electron-electron interactions |38l Il2])- In the perturbative limit the self- 
energies can be calculated up to order A^ using the Keldysh formalism. The second order 
expressions are 



E(f (c.) = J ^[G^r (c - z.)Z^(°)+-(z.) + ^^"+(0; - u)D('^^{u)] (35) 

where D^^\uj) is the unperturbed phonon propagator and G'^q{u}) are the electronic propa- 
gators of the quantum level for the nonequilibrium effective one electron problem. 



From Eqs. (35) it is straightforward to verify that in the limit F — t- 0, T,^fj^{u) tends to 



an expression formally identical to that of Eq. (33) in the equilibrium situation. Therefore 



the ansatz of Eq. (34) will recover automatically i) the atomic limit and ii) the results of 
nonequilibrium perturbation theory in the limit A/F — > 0. There still remains the problem 
of finding an analogous interpolative ansatz for the Keldysh self-energies and S 
[55| I12H11] . This is not as straightforward as in the retarded case because these self-energies 
are not well defined in the atomic limit. An appropriate ansatz can however be obtained by 
requiring that T,^~ and S ^ satisfy the Keldysh relation [381111]: 

SoV(c^) - Soo+(u;) = S5o(^) - SSo(^) = 2^/mSSo(a;) (36) 

and that the results of second order perturbation theory will be recovered in the limit 
A/F — )■ 0. This conditions are fulfilled by the following ansatz: 

JmS^o^ (w) 

with an analogous expression for T,Qf^{u)). In addition to the above requirements this ex- 
pression recovers the equilibrium limit: 

S+-(u;) = 2z/mGSo(c^)/H (38) 
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where /(w) is the equihbrium Fermi distribution function of the electrodes. An analogous 



ansatz to the one of Eq. (37) was used in [3H1I12] for the case of electron-electron interactions. 

Finally we will comment on the self-consistency procedure. We impose consistency in 
the dot charge between the one-electron and the interacting problem. This is achieved by 
introducing an effective dot level position in the one-electron Hamiltonian. As we mentioned 
at the end of section III, current conservation is not necessarily fulfilled for an approximate 
solution. In particular this is the case for the ISA. The self-consistent procedure can be 
nevertheless generalized by requiring both charge and current consistency between the one- 
electron and interacting cases |3Sj. A natural choice is to introduce effective one-electron 
couplings of the dot level with the electrodes T^^eff, ^R,eff which are fixed from the require- 
ment of current consistency. When imposing current consistency, the current is calculated 
by means of Eqs.([8]); otherwise Eq.([9]) is used. In this work and for the range of parameters 
considered, the requirement of current consistency does not alter in a significant way the 
results for current and conductance but the agreement with the EOM results somewhat 
improves when imposing it. 



V. RESULTS AND DISCUSSION 

All the calculations in this work have been performed in the limit of zero temperature so 
only phonon emission is possible. Currents are plotted in units of e/h and conductances are 
plotted in units of e'^/h . 

It will be useful for the discussion of our results to have a scheme of the inelastic processes 
we describe in this work. Fig. la sketches a process in which an electron from the left 
electrode tunnels to the localized level where it emits n phonons. Energy conservation 
requires /iL — cq = uuq (or — fi^ = tiuq). The threshold for emission of n phonons is 
depicted in Fig. lb. It occurs when an electron from the left electrode jumps into the right 
electrode through the level and therefore requires V = fiL — fiR = nojQ. Both processes show 
up in the conductance of the system with characteristic signatures that we analyze in this 
section. 

We start this section by discussing a situation in which the bias potential V is applied 
symmetrically between the electrodes so that fin = —fiL = —eV/2. The Fermi energy of the 
leads in equilibrium is taken as our zero of energy. The symmetry of the problem makes 
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FIG. 1. (Color online) Sketch of the inelastic processes we analyze in this work. In (a) an electron 
from the left electrode tunnels to level, where it emits n phonons, and passes to the right electrode. 
The onset for emission on n phonons is illustrated in (b), where one electron at the left chemical 
potential tunnels to the level, emits n phonons and continues at the right chemical potential. Blue 
thin arrows represent tunneling events and thick arrows phonon emission. 

the I — V curves and the conductance to be identical for negative and positive values of eg. 
Also I{—V) = —I{V). Then we show results only for positive values of eg and V. 

Fig. 2 shows the current (upper panels) and the conductance (lower panels), as a function 
of V for A = O.Suo, F = O.Iwq and for three values of the gate potential corresponding to 
eo = 0.0, 0.3 andO.ScjQ. The current and the conductance are separated into their elastic 
and inelastic contributions, showing clearly that, for this small value of X/ujq, the current 
is predominantly elastic. The inelastic current has a threshold at the onset for inelastic 
processes, V = for phonon emission which shows up as a step in the conductance. Even 
though this step is tiny on the scale of this figure in cases (a) and (b) because of the small 
value of X/uq used here, it is a feature that we will discuss extensively in the context of 
fig.4. The conductance also shows different lorentzian-like peaks at V = 2|eo ± ^i^ol, with 
n a positive integer. These peaks are the signature of the inelastic processes described in 
Fig.la. For the case eo = of Fig.2a they appear at = 0,2wo, 4cJo ••• while for cq > 
each peak is split into two which, according to the energy conservation requirements stated 
above, appear at \^ = 2|eo ± nuo\. The peaks of the conductance correspond to the steps in 
the I — V curves, their width being proportional to F. 

Fig. 3 is as Fig. 2 but we have increased the value of the electron-phonon interaction to 
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FIG. 2. (Color online) Currents (upper panels) and conductances (lower panels) as a function of 
V (in units of loq) for a symmetrically applied bias voltage and X/ojo = 0.3 . Red lines: elastic 
components, green lines: inelastic components, black lines: total values. Continuous lines: EOM, 
dotted lines: ISA. (a) eo = 0, (b) eo = O.Suq and (c): eo = O.Swo 

A = O.Two (while keeping the same values of eo). For this value of A we are far from the 
perturbative regime and the steps in the I — V curve for n = 2 are clearly visible. Notice how 
the contribution of the inelastic processes to the total current and the conductance increases 
quickly with the applied bias, overcoming the contribution of the elastic processes, as we 
move away from the electron-hole symmetric case eo = 0. This behavior, in which the current 
versus voltage curves tend to adopt a staircase form with steps located at = 2|eo ± nuo\, 
is enhanced as X/uq gets larger than 1. The height of the steps in the current gives the 
probability of emitting n phonons and follows very approximately the Poisson distribution, 
e~^^, with g = (^)^- This behavior is qualitatively similar to what was obtained in Ref.|llj 
using a semiclassical master equations approach. The staircase behavior of conductance with 
applied bias due to phonon emission has been experimentaly found in Ref|5j. The main 
peak of the conductance, obtained at low voltages, is almost identical to the main resonance 
appearing in the equilibrium density of states, showing the polaronic reduction of the level 
width [8]. However, the phonon side-bands show specific features associated to inelastic 
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FIG. 3. (Color online) Currents (upper panels) and conductances (lower panels) as a function of 
V (in units of loq) for a symmetrically applied bias voltage and A/wq = 0.7 . Red lines: elastic 
components, green lines: inelastic components, black lines: total values. Continuous lines: EOM, 
dotted lines: ISA. (a) eo = 0, (b) eo = O.Sojq and (c): eo = O.Bloq 



transport, which we will analyze next. The total conductance shows steps at = tiuq. We 
should mention that not only the inelastic component exhibits this feature but the elastic 
component as well because of the change in the retarded self-energy due to the appearance 
of new inelastic processes. 

In Figs. 2 and 3 we compare the results from both theoretical approaches, EOM and 
ISA. The remarkable agreement found gives confidence in the interpolative scheme and also 
in the EOM method to the order 0{V^^) for values of X/uq up to 1. At this point we 
should comment that the EOM method up to the order Oiy^^^) starts to show numerical 
instabilities for higher values of A associated with the increasing number of phonons that have 



to be included in the solution of Eqs. (16) and (26) and with the corresponding logarithmic 
singularities in F'-'^-*''-^^ (Eq.(17)). This problem was already found in equilibrium and it is 
cured by the renormalization of these singularities that appears when the method is carried 
to the order Oiy^^^). However, the extension of the procedure to situations out of equilibrium 
is not straightforward and will be deferred to further work. 
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FIG. 4. (Color online) The conductance as a function of V (in units of ujq) for a symmetrically 
applied voltage and X/ujq = 0.5. Panels (a), (b) and (c) show the regions near V = ujq, V = 2u;o 
and V = 3uJo respectively. The results of the EOM method are shown for: eo = (black lines), 
eo = O.lwo (red lines), eo = 0.2a;o (green lines), eo = O.Swo (blue lines) eo = OAcoq (magenta lines), 
eo = 0.5a;o (orange lines). 

As mentioned in the Introduction, the issue of whether the steps in the total conductance 
a.t V = are upwards or downwards has raised a great interest both theoretically and 
experimentally. Both our formalisms recover the results already obtained in the weak cou- 
pling regime and in the following we concentrate in the regime of strong coupling, X,uo > T, 
where we find jumps of the conductance aX V = nuo for any n. 

Fig. 4 shows the conductance as a function of the applied bias voltage in the regions 
near: (a) ojq, (b) 2uo and (c) 3wo, for F = O.Icjq, A = O.Swq and several values of the gate 
voltage corresponding to eo = 0,0.1,0.2,0.3,0.4 and O.Swo . For the sake of clarity, only 
the results of the calculations using the EOM method are shown. Note in Figs. 4 (a) and 
(c) that the step in the conductance is always upwards except for eo = O.Scuo, where it is 
downwards and the conductance is at a relative maximum. The same happens in Fig. 4(b), 
with the conductance jumping downwards only for cq = 0, for which value the conductance 
has a relative maximum at \^ = 2ujq . These results can be understood in terms of the 
interference between the step-like processes at = nuo and the lorentzian-like peaks at 
V = 2|eo in'cuol- When both conditions do not coincide, the inelastic conductance increases 
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at V = nojQ and dominates the elastic decrease which is very small there. Consequently, 
the conductance step is upwards. However, if 2|eo ± n'ujQ\ = uujq (within an accuracy of 
±r), we always find a downward decrease of the total conductance steps. The origin of 
this behavior is different for n = 1 than for the rest of the cases. The value V = uq is the 
absolute onset for inelastic processes and, consequently, the inelastic conductance increases 
there. This increase is compensated by a stronger decrease of the elastic conductance in a 
way similar to the one analyzed theoretically in the perturbative regime F >> A, Uq [261 - I28] . 
However, for n > 1 we find the inelastic conductance decreasing at V = huq while the elastic 
one increases there. The appearance of a new inelastic channel emitting n phonons makes 
the intensity of the previously existing ones to decrease abruptly. Thus we attribute the 
different behaviors of the elastic/inelastic components of the conductance to interferences 
between the inelastic processes of Figs, la and b, which can occur for n = 2, 3.... The total 
conductance always shows a downward step whenever the value V = huq is at a relative 
maximum. In any other case, the conductance jumps up at = tiuq. This seems to be a 
very general behavior, valid in both the strong and weak coupling regimes in A/F, at least 
in cases of symmetric coupling between the localized level and the electrodes. It is seen for 
any value of A not only for = 1, as the perturbation theory predicts, but for any value of 
n. 

We have already pointed out the good agreement obtained by our two theoretical ap- 
proaches in the case of a symmetrically applied bias. That this agreement is not fortuitous 
is proved by comparing the results in a different situation, in which the bias is applied 
asymmetrically, with fipt = ^ and fi^ = V- This is done in Figs. 5 and 6, where we show the 
current and the conductance for A = O.Suq and A = O.Twq respectively, for several values of 
eo. For simphcity, we have chosen F^ = F^j = F/2 with F = O.Icjq. Only positive values 
of eo are shown because /(— eo,^) = — -^(^o, — ^) and G{—eQ,V) = G{eo,—V). At variance 
from Figs. 2 and 3, the maximum of the conductance is very close to 1. The larger deviations 
from perfect conductance are obtained in Fig.6, for large eo which means that we are far 
from equilibrium. The fact that the EOM results are higher than the interpolative results 
at the maximum is the consequence of the numerical inaccuracies commented above. As 
in Figs. 2 and 3, the current increases in a step- like way. Correspondingly, the conductance 
presents lorentzian-like phonon side-bands associated with the inelastic process occurring at 
V = io ± nuo and jumps at V = muo, with a strong change in line shape under conditions 
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FIG. 5. (Color online) Currents (upper panel) and conductances (lower panel) as a function of 
V = fiL (in units of wq) for X/ujq = 0.3 and: eo = (black lines), eo = O.Iujq (red lines), eo = 0.2u;o 
(green lines), eo = O.Scjq (blue lines) and eo = 0.5ujo (magenta lines). Continuous lines: EOM, 
dotted lines: ISA. 




when they can both occur and interfere. Therefore, this is a robust behavior obtained by 
both theoretical approaches under different values of the parameters defining the problem. 
The asymmetry of the conductance for positive and negative values of V is a consequence 
of the very asymmetric behavior of the level occupancy < nQ{V) > when one of the elec- 
trodes do not change its chemical potential. This can be qualitatively understood from the 



atomic Green function, Eq. (30), where one can readily see that, for positive values of eo 
and Uo > eo, phonon emission with V > {V < 0) should be proportional to 1— < no > 
(< no >). Also, the asymmetry of the conductance follows the shape of the nonequilib- 
rium density of states (not shown) with V > (V < 0) mapping out its empty (occupied) 
portions. 



VI. CONCLUSIONS 



In this work, we present a combined theoretical approach to analyze the nonequilibrium 
transport properties of nanoscale systems exhibiting strong electron-phonon interactions and 
coupled to metallic electrodes. We describe the system by the spinless Anderson-Holstein 
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FIG. 6. (Color online) Currents (upper panel) and conductances (lower panel) as a function of 
V = fiL (in units of wq) for X/ujq = 0.7 and: eo = (black lines), eo = O.Iujq (red lines), eo = O.Swq 
(green lines) Continuous lines: EOM, dotted lines: ISA. 



Hamiltonian and use a Keldysh Green function formalism to generalize an equation of motion 
method and an interpolative self-energy approximation previously developed in equilibrium. 
These two approaches recover the results obtained formerly in the weak coupling regime 
A/r << 1 and this article concentrates in the strong coupling regime A,a;o > T. Using both 
techniques, we analyze the specific features of inelastic transport appearing in the intensity 
versus voltage curves and in the conductance. Excellent overall agreement between both 
approaches is found in a wide range of parameters. We obtain a step-like increase of the 
current with the applied voltage at V = eo± nuo with the corresponding phonon sidebands 
of the conductance, a behavior which gets more pronounced as A/cuq increases. We also find 
steps in the conductance at = nuo for any value of n. These are generally upwards, except 
when the value V = nuo occurs at a relative maximum of the conductance in which case it 
is downwards. This seems to be a very general behavior, valid in both the strong and weak 
coupling regimes in A/F, at least in cases of symmetric coupling between the localized level 
and the electrodes. 
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